%Integration
Now that the dynamics of the system are defined we want to be able to integrate over the time domain to actually simulate the deformation. This can be done by either using explicit or implicit integration. Explicit integration has the upside of being easy to implement, however it suffers from being only conditionally stable and the errors accumulate. Implicit integration on the other hand can be slightly harder to work with but is instead unconditionally stable thus recommended for simulating this kind of system. 

Since the dynamic system was described in the previous section, the equation of motion \ref{eqn:fullDiscPDE} can now be applied to the Implicit Euler form as follows \cite{muller_ivm}

\begin{equation}\label{eqn:impEuler}
		(\mathbf{M} + \Delta t  \mathbf{C} + \Delta t^2 \mathbf{K})v^{t+1} = \mathbf{M} v - \Delta t (\mathbf{K}(\mathbf{x} - \mathbf{x}_0) +\mathbf{f}_{ext})
\end{equation}

Where the velocity for the next time step $v^{t+1}$ is what needs to be solved. It can be noticed that equation \ref{eqn:impEuler} is in the form of the discrete Poisson equation $\mathbf{Ax=b}$ which means that the inverse of $\mathbf{A}$ has to be computed. Since this is a rather computationally heavy procedure for especially large matrices, it can be a good idea to consider the use of an algorithm that approximates the result of the velocity, like for example the Conjugate Gradient method.

The Conjugate Gradient method is an iterative numerical algorithm that requires the matrix $\mathbf{A}$ in the discrete Poisson equation to be both sparse and symmetrical. Basically what it does is that it takes the equation as input and solves $\mathbf{x}$, thus making it very applicable for this kind of problem.